Finite difference method for transport properties of massless Dirac fermions 
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We adapt a finite difference method of solution of the two-dimensional massless Dirac equation, 
developed in the context of lattice gauge theory, to the calculation of electrical conduction in a 
graphene sheet or on the surface of a topological insulator. The discretized Dirac equation retains a 
single Dirac point (no "fermion doubling"), avoids intervalley scattering as well as trigonal warping, 
and preserves the single- valley time reversal symmetry {— symplectic symmetry) at all length scales 
and energies — at the expense of a nonlocal finite difference approximation of the differential 
operator. We demonstrate the symplectic symmetry by calculating the scaling of the conductivity 
with sample size, obtaining the logarithmic increase due to antilocalization. We also calculate 
the sample-to-sample conductance fluctuations as well as the shot noise power, and compare with 
analytical predictions. 

PACS numbers; 71.10.Fd, 73.20.Fz,73.23.-b 



CD 



o 



> 

OO 

d 

00 

o 



X 



I. INTRODUCTION 

The discovery of graphene [J] has created a need for 
efficient numerical methods to calculate transport prop- 
erties of massless Dirac fermions. The two-dimensional 
massless Dirac equation (or Weyl equation) that governs 
the low-energy and long- wave length dynamics of conduc- 
tion electrons in graphene has a time reversal symmetry 
called symplectic — which is special because it squares to 
— 1. (The usual time reversal symmetry, which squares 
to is called orthogonal.) The symplectic symmetry 
is at the origin of some of the unusual transport proper- 
ties of graphene Jl, H, 0, H, @] , including the absence of 
back scattering [J, weak antilocalization enhanced 
conductance fluctuations [1, [lOj , and absence of a metal- 
insulator transition pdl . [T3 |. 

Numerical methods of solution can be divided into two 
classes, depending on whether they break or preserve the 
symplectic symmetry. 

The tight-binding model of graphene, with nearest 
neighbor hopping on a honeycomb lattice, breaks the 
symplectic symmetry by the two mechanisms of interval- 
ley scattering [8| and trigonal warping [13i] . Intervalley 
scattering couples the two flavors of Dirac fermions, cor- 
responding to the two different valleys (at opposite cor- 
ners of the Brillouin zone) in the graphene band struc- 
ture, thereby changing the symmetry class from sym- 
plectic to orthogonal. Trigonal warping is a triangular 
distortion of the conical band structure that breaks the 
momentum inversion symmetry {-\-p —>■ — p), thereby ef- 
fectively breaking time reversal symmetry in a single val- 
ley and changing the symmetry class from symplectic to 
unitary. 

Breaking of the symplectic symmetry eliminates both 
weak antilocalization as well as the enhancement of the 
conductance fluctuations, and drives the system to an 
insulator with increasing size or disorder n^. ITU . As 
observed in computer simulations [l^, llSj. the break- 
ing of the symplectic symmetry can be pushed to larger 



system sizes and larger disorder strengths by reducing 
the lattice constant (at fixed correlation length and fixed 
amplitude of the disorder potential) — but this severely 
limits the computational efficiency. 

The Chalker-Coddington network model ^ HO, HI, 
applied to graphene in Ref. [22], has a single flavor of 
Dirac fermions, so there is no intervalley scattering — 
but it still belongs to the same class of methods that 
break the symplectic symmetry of the massless Dirac 
equation. (The symplectic symmetry is broken on short 
length scales by the Aharonov-Bohm phases that appear 
in the mapping of the Dirac equation onto the network 
model.) 

Both the network model and the tight-binding model 
are real space regularizations of the Dirac equation, with 
a smallest length scale (the lattice constant) to cut off the 
unbounded spectrum at large positive and large negative 
energies. There exists at present just one method to cal- 
culate transport properties numerically while preserving 
the symplectic symmetry, developed independently (and 
implemented differently) in Refs. [ll| and [l^l- That 
method (used also in Refs. [l^, H^l) is based on a mo- 
mentum space regularization, with a cutoff of the Fourier 
transformed Dirac equation at some large value of mo- 
mentum. 

It is the purpose of the present paper to develop and 
implement an alternative method of solution of the Dirac 
equation, that shares with the tight-binding and network 
models the convenience of a formulation in real space 
rather than momentum space, but without breaking the 
symplectic symmetry. 

A celebrated no-go theorem [25] in lattice gauge theory 
forbids any regularization of the Dirac equation with local 
couplings from preserving symplectic symmetry. (The 
problematic role of intervalley scattering appears in that 
context as the fermion doubling problem.) Several nonlo- 
cal finite difference methods have been proposed to work 
around the no-go theorem and we will adapt one of these 
(developed by Stacey [2^ and by Bender, Milton, and 
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Sharp [131) to the study of transport properties. 

The adaptation amounts to 1) the inclusion of a spa- 
tially dependent electrostatic potential (which breaks 
the chiral symmetry that played a central role in Rcfs. 
[H, [13|)> ^^'^ 2) a proper discretization of the current 
operator (such that the total current through any cross 
section is conserved). We implement the finite difference 
method to solve the scattering problem of Dirac fermions 
in a disordered potential landscape connected to ballistic 
leads, and compare our numerical results for the scaling 
and statistics of conductance and shot noise power with 
analytical theories [§, [13, HI] . 

Our numerical method is relevant for electrical con- 
duction in graphene under the assumption that the im- 
purity potential in the carbon monolayer is long-ranged 
(so that intervalley scattering is suppressed) and weak 
(so that trigonal warping can be neglected). Massless 
Dirac fermions are also expected to govern the electrical 
conduction along the surface of a three-dimensional topo- 
logical insulator [sO, HH (recently realized in BiSb 
[32|). In that case the symplectic symmetry is preserved 
even for short-range scatterers, and our numerical results 
should be applicable more generally. 

II. FINITE DIFFERENCE REPRESENTATION 
OF THE TRANSFER MATRIX 

A. Dirac equation 

We consider the two-dimensional massless Dirac equa- 
tion, 

H^ = E^, H = ~ihv{(TA + (Tydy) + U{r), (2.1) 

where v and E are the velocity and energy of the Dirac 
fermions, U{x,y) is the electrostatic potential landscape, 
and ^(x, y) is the two-component (spinor) wave function. 
The two spinor components of ^' refer to the two atoms 
in the unit cell in the application to graphene, or to the 
two spin degrees of freedom in the application to the 
surface of a topological insulator. We note the symplectic 
symmetry of the massless Dirac Hamiltonian, 

H ^ SHS^^ = <TyH*ay. (2.2) 

The time reversal symmetry operator S — iUyC (with C 
the operator of complex conjugation) squares to —1. The 
chiral symmetry UzHcTz = —H is broken by a nonzero U, 
so it will play no role in what follows. For later use we 
also note the current operator 

{jx,jy) = v{a^,ay). (2.3) 

We consider a strip geometry of length L along the lon- 
gitudinal X direction and width W along the transversal 
y direction. For the discretization we use a square lat- 
tice, Xm = mA, yn = nA, with indices m = 1,2, ... AI 
(M = L/A), n = 1,2,...N {N = W/A). In the ap- 
plications we will consider large aspect ratios W/L ^ 1, 



for which the precise choice of boundary conditions in the 
transverse direction does not matter. We choose periodic 
boundary conditions, yN+i = yi, since they preserve the 
symplectic symmetry. The values '^m,n — '^{xrmyn) of 
the wave function at a lattice point are collected into a set 
of iV-component vectors = (*m,i, ^'m,2, ■ • • ^'m.w)'^, 
one for each m = 1, 2, . . . M. 

The N X N transfer matrix Aim is defined by 

*m+l = Xm*m- (2.4) 

The symplectic symmetry (|2.2p of the Hamiltonian re- 
quires that 'J and a-y'i'* are both solutions at the same 
energy E, so they should both satisfy Eq. (|2.4p . The 
corresponding condition on the transfer matrix is 

M,n = cryMmay. (2.5) 

The transfer matrix should conserve the total current 
through any cross section of the strip. In terms of the 
(still to be determined) discretized current operator J^, 
this condition reads {^m+i\Jx\'^m+i) = (*,„| Jxl*™), 
which then corresponds to the following condition on the 
transfer matrix: 

Mlj^Mrn = Jx. (2.6) 

Our problem is to discretize the differential operators 
in the Dirac equation (|2.ip . as well as the current op- 
erator (|2.3p . in such a way that the resulting transfer 
matrix describes a single flavor of Dirac fermions and 
without violating the two conditions (|2.5p . (|2.6p of sym- 
plectic symmetry and current conservation. 



B. Discretization 

A local replacement of the differential operators dx , dy 
by finite differences either violates the Hermiticity of H 
(thus violating conservation of current) or breaks the 
symplectic symmetry (by the mechanism of fermion dou- 
bling) . A nonlocal finite difference method that preserves 
the Hermiticity and symplectic symmetry of H was devel- 
oped by Stacey [2^ and by Bender, Milton, and Sharp 
j27| . These authors considered the case U — 0, when 
both symplectic and chiral symmetry are present. We 
extend their method to a spatially dependent U (thereby 
breaking the chiral symmetry) , and obtain the discretized 
transfer matrix and current operator. 

Since the transfer matrix relates ^(x,?/) at two differ- 
ent values of x, it is convenient to isolate the derivative 
with respect to x from the Dirac equation (|2.ip . Multi- 
plication of both sides by (i/hv)ax gives 

dx"^ = i-ia,dy - laxV)^, (2.7) 

with the definition V — {U — E)/hv. We can now make 
contact with the discretization in Refs. [1^ [131 of the 
Dirac equation in one space and one time dimension, with 
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FIG. 1: Square lattice (filled circles) on which the wave func- 
tion is discretized as ^m.n- The finite differences are eval- 
uated at the displaced points indicated by crosses. The Dirac 
equation (|2.7p is applied at the empty circles, by taking the 
mean of the contributions from the two adjacent crosses. The 
resulting finite difference equation defines a transfer matrix 
in the x direction that conserves current and preserves the 
symplectic symmetry. 



that J and V*^™^ are real symmetric matrices, while K, 
is real antisymmetric. Furthermore J and K. commute, 
but neither matrix commutes with V'™-'. 

For later use, we note that J has eigenvalues 



ji =2cos^{Trl/N), ; = 1,2,...7V, 



(2.14) 



corresponding to the eigenvectors V''^'-' with elements 

i^ii'^ = N-^/^exp{27riln/N). (2.15) 

The eigenvalues of IC are 

Ki=ism{2nl/N), l = l,2,...N, (2.16) 

for the same eigenvectors tp'^^^ From Eq. (|2.14p we see 
that for N even there is a zero eigenvalue oi (at I = 
N/2). To avoid the complications from a noninvertible 
J, we restrict ourselves to N odd (when all eigenvalues 
of J7 are nonzero). 

C. Transfer matrix 



X playing the role of (imaginary) time and y being the 
spatial dimension. 

The key step by which Refs. [1^, [131 avoid fermion 
doubling is the evaluation of the finite differences on a 
lattice that is displaced symmetrically from the original 
lattice. The displaced lattice points [xm + A/2,y„ -|- 
A/2) are indicated by crosses in Fig.[TJ On the displaced 
lattice, the differential operators are discretized by 

(2.8) 
(2.9) 

and the potential term is replaced by 

(2.10) 

with Vm.n — ^(a;m + A/2, y„-|- A/2). The Dirac equation 
()2.7|) is applied at the points {xm + A/2, ?/„) (empty cir- 
cles in Fig.[T]) by averaging the terms at the two adjacent 
points {x,n + A/2, y„ ± A/2). 

The resulting finite difference equation can be written 
in a compact form with the help of the N x N tridiagonal 
matrices J, K., V*-™-*, defined by the following nonzero 
elements: 

\Jn,n — 1; yJn.n-\-\ — Un^n — X — 2' (^•-^-^) 

■>;(m) _ in/ -L V \ -i;'-™) — It/ 

^ntl = hVm.n-l- (2.13) 

In accordance with the periodic boundary conditions, the 
indices n ± 1 should be evaluated modulo N. Notice 



The discretized Dirac equation is expressed in terms of 
the matrices (ICT1) - (PT^ by 

(*,„ + *„,+i). (2.17) 

Rearranging Eq. (|2.17p we arrive at Eq. (|2.4p with the 
transfer matrix 

M,n= (j + M + izAa,V(™))"' 

( J - ia.lC - izAa^V^"') . (2.18) 

Since we take N odd, so that J' is invertible, we may 
equivalently write Eq. (|2.18p in the more compact form 

(2.19) 

As announced, the transfer matrix is nonlocal (in the 
sense that multiplication of '5'™ by couples all trans- 
verse coordinates). 

One can readily check that the condition (|2.5p of sym- 
plectic symmetry is fuUfilled. In App. lAl we demonstrate 
that the condition (|2.6p of current conservation holds if 
we define the discretized current operator in terms of 
the symmetric matrix ^7, 

Jx = ^va.xJ. (2.20) 

The absence of fermion doubling is checked in Sec. IIV Al 
The transfer matrix Ai through the entire strip (from 
a; = to a; = L) is the product of the one-step transfer 
matrices Mm, 

M 

M=l[ Mm, (2.21) 

m— 1 
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ordered such that Aim+i is to the left of Mm- The prop- 
erties of symplectic symmetry and current conservation 
are preserved upon matrix multiphcation. 

D. Numerical stability 

The repeated multiplication (|2.21|1 of the onc-stcp 
transfer matrix to arrive at the transfer matrix of the 
entire strip is unstable because it produces both exponen- 
tially growing and exponentially decaying eigenvalues, 
and the limited numerical accuracy prevents one from 
retaining both sets of eigenvalues. We resolve this obsta- 
cle, following Refs. [iTI, 123, [s^l , by converting the transfer 
matrix into a unitary matrix, which has only eigenvalues 
of unit absolute value. The formulas that accomplish this 
transformation are given in App. [Bl 

III. FROM TRANSFER MATRIX TO 
SCATTERING MATRIX AND CONDUCTANCE 

A. General formulation 

The scattering matrix is obtained from the transfer 
matrix by connecting the two ends of the strip aX x — 
and X — L to semi-infinite ballistic leads. The N trans- 
verse modes in the leads (calculated in Sec. lIVp . consist 
of A^o propagating modes (fy^ (labeled -I- for right-moving 
and — for left-moving), and N — Nq evanescent modes 
xf^ (decaying for x ±oo). The propagating modes are 
normalized such that each carries unit current. 

Consider an incoming wave in mode Iq from the left. 
At a; = 0, the sum of incoming, reflected, and evanescent 
waves is given by 

<* = -^/l; + E^Mo0r + E"'.'oXz-, (3.1) 
I I 

while the sum of transmitted and evanescent waves at 
a; = L is given by 

'i>;^ = E*'.'o0z++E"Ux^ (3.2) 

The Nq X Nq reflection matrix r and transmission ma- 
trix t are obtained by equating 

<i>[f* = A^$|f, (3.3) 

eliminating the coefficients a, a', and repeating for each 
of the Nq propagating modes incident from the left. 
Starting from a mode incident from the right, we sim- 
ilarly obtain the reflection matrices r' and t' , which to- 
gether with r and t form a 2A'o x 2Ao unitary scattering 
matrix, 

5= (;:;). (3.4) 



As a consequence of unitarity, the matrix products tt^ 
and t'i'^ have the same set of eigenvalues Ti , T2 , . . . T^o , 
called transmission eigenvalues. 

The number Nq of propagating modes in the leads is 
an odd integer, because of our choice of periodic bound- 
ary conditions. The symplectic symmetry condition (j2.5p 
then implies that the transmission eigenvalues r„ consist 
of one unit eigenvalue and {Nq — l)/2 degenerate pairs 
(Kramers degeneracy 0|). 

The conductance G follows from the transmission 
eigenvalues via the Landauer formula, 

G-GoE^"- (3.5) 

n 

The conductance quantum Gq = 4e^/ft, in the application 
to graphene (which has both spin and valley degenera- 
cies), while Go = e^/hin the application to the surface of 
a topological insulator. The Kramers degeneracy, which 
is present in both applications, is accounted for in the 
sum over the transmission eigenvalues. 

B. Infinite wave vector limit 

Following Ref. [s^, we model metal contacts by leads 
with an infinitely large Fermi wave vector. In the infinite 
wave vector limit all modes in the leads are propagating, 
so Nq = N and the scattering matrix has dimension 2 A^ x 
2N. The states <j)f {I = 1,2,...A^) in this limit are 
simply the 2N eigenstates of the current operator J^, 
normalized such that each carries the same current. In 
terms of the eigenvalues and eigenvectors (12.141) . (|2.15p 
of J' we have 

'^^-jT^/^(±\)^«. (3.6) 

Instead of the general Eqs. (|3.ip and (|3.2p we now have 
the simpler equations 

= '/'^+E-Mo'^r, 'i>;r = E^Mo'^z+. (3.7) 

To obtain from Eq. p.3p a closed-form expression for S 
in terms of M, we first perform the similarity transfor- 
mation 

M = nMTZ-\ 'R = (thJ^'^, (3.8) 
where an is the Hadamard matrix, 

a^=2-V2Q l^^^.^i. (3.9) 

The notation ohJ^^"^ signifies a direct product, where 
an acts on the spinor degrees of freedom s = ± and J^^"^ 
acts on the lattice degrees of freedom n = 1, 2, . . . A^. No- 
tice that the matrix 7?. is Hermitian [since J is Hermitian 
with exclusively positive eigenvalues, see Eq. (|2.14p ]. 
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We separate the spinor degrees of freedom of into 
four N X N blocks, 



M 



M- 



M 
M 



/ ' 



(3.10) 



such that Mns,ms' = Mnm- The matrix M. has a corre- 
sponding decomposition into submatrices J^/i"^ . As one 
can verify by substitution into Eq. (|3.7p and comparison 
with Eq. (jS.Sp . the submatrices M.'^'^ are related to the 
transmission and reflection matrices by 




r = -[M 



t = 

t' = i^M- 
r' = 



-1 , 



M- 



(3.11a) 
(3.11b) 
(3.11c) 
(3. lid) 



Similar formulas were derived in Ref. [22| , but there the 
transformation from to S" involved only a Hadamard 
matrix and no matrix J ^ because of the different current 
operator in that model. 



IV. BALLISTIC TRANSPORT 

For a constant U we have ballistic transport through 
the strip of length L and width W . In this section we 
check that we recover the known results [s^ for ballis- 
tic transport of Dirac fermions from the discretized trans- 
fer matrix. 



A. Dispersion relation 



For U — Uq — constant the matrix (|2.13p of discretized 
potentials is given by V'-™^ = —{e/A)J, with e — [E — 
UQ)A/hv the dimensionless energy (measured relative to 
the Dirac point at energy Uq). Substitution into Eq. 
(|2.19p gives the m-independcnt ballistic transfer matrix 



Xball = 



(4.1) 



This is the one-step transfer matrix. The transfer matrix 
through the entire strip, in this ballistic case, is simply 

M = (A^ball)*'. 

In accordance with Eqs. (|2.14p - (|2.16p . the matrix 
J^^JC can be diagonalized (for N odd) by 

J^^1C = A„„' = i tan(7rn/7V)(5„„/ , (4.2a) 

= iV-i/2 exp(27rmn7^). (4.2b) 



FIG. 2: Dispersion relation (|4.5|) of the discretized Dirac equa- 
tion, plotted in the first Brillouin zone for two transverse 
modes (I — N, solid curve; I ~ A^/4, dotted curve). The dis- 
persion relation approaches that of the Dirac equation near 
the point (fc, e) = (0, 0), and avoids fermion doubling at other 
points in the Brillouin zone. 



The Fourier transformed transfer matrix J-A4hediJ^^ is 
diagonal in the mode index Z = l,2,...Af. A2x2 matrix 
structure mi in the spin index remains, given by 



mi 



1 + i{e/2)a,: + ta.n{-Kl/N)(j-, 
1 - i{e/2)a^ - tan(7r?/iV)cr^ 

The eigenvalues and eigenvectors of mi are 



(4.3) 



e/2 

iia,-a{-Kl/N) ±tan(fc/2) ) ' 

(4-4) 

with the dimensionless momentum k given as a function 
of e and I by the dispersion relation 



tan2(fc/2) +tan2(7rVA^) = (e/2)" 



(4.5) 



In Fig. [5] we have plotted the dispersion relation (|4.5p 
for two different modes in the first Brillouin zone — tt < 
k < TT. For each mode index I, there is one wave that 
propagates to positive x (on the branch with de/dk > 
0) and one wave that propagates to negative x (on the 
branch with de/dk < 0]. 

As anticipated [l^, [13] , the discretization of the Dirac 
equation on the displaced lattice (crosses in Fig. [T|) has 
avoided the spurious doubling of the fermion degrees of 
freedom that would have happened if the finite differ- 
ences would have been calculated on the original lattice 
(solid dots in Fig. [1]). In the low-energy and long- wave- 
length limit fc,e — > 0, the conical dispersion relation 
(vpx)'^ + {vPyY = (E — Uq)'^ of the Dirac equation (12. ip 
is recovered. The longitudinal momentum is px = hk/ A, 
while the transverse momentum is py = {2t:1i/W)1 if 
l/N ^Qovpy^ -{2TTh/W){N - I) if l/N 1. 



B. Evanescent modes 

For |e| < 2tan(7r/A^), hence for \E - Uq\ < 2nnv/W, 
only the mode with index I = N is propagating. The 
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FIG. 3: Relation between the energy e and the imaginary part 
K of the wave number of evanescent modes, calculated from 
Eq. (14. 6p for five different values of the mode index [parame- 
terized by ^ = tan(7rZ/A'^)]. The real part of the wave number 
equals on the solid contours (corresponding to «:+), while it 
equals tt on the dashed contours (corresponding to Only 
the K+ evanescent modes have a correspondence to the Dirac 
equation in the limit e — > 0. The K- evanescent modes that 
appear for |5| > 1 are artefacts of the discretization for large 
transverse momenta. 



other — 1 modes are evanescent, that is to say, their 
wave number k has a nonzero imaginary part k. There 
are two classes of evanescent modes, one class with a 
purely imaginary wave number k = ik^, and another 
class with a complex wave number fc = vr + The 
relation between k± and e, following from Eq. (|4.5p . is 

tanh (k+/2) = tan2(7r//A^) - (e/2)^ (4.6a) 
cotanh (k_/2) tan2(7r//A^) - ie/2f. (4.6b) 

In Fig. [3] we have plotted Eq. (|4.6p for different mode 
indices, parameterized by ^ = tan(7rZ/A^). The evanes- 
cent modes in the Dirac equation correspond to fc = iK+ 
in the limit e ^ (solid contours in Fig.[3|). The second 
"spurious" class of evanescent modes, with k — n + in^ 
(dashed contours), is an artefact of the discretization 
that appears for large transverse momenta (|^| > 1, or 
iV/4 < / < 3iV/4). 

To minimize the effect of the spurious evanescent 
modes we insert a pair of filters of length Lq between the 
strip of length L and the leads with infinitely large Fermi 
wave vector. By choosing a large but finite Fermi wave 
vector in the filters, they remove the spurious evanescent 
modes of large transverse momenta which are excited by 
the infinite Fermi wave vector in the leads. 

The geometry is sketched in Fig. 2] In the filters we 
choose U = E — 2fiv/lS. (so e = 2 in the filters). 
Since > tt/A^ for the spurious evanescent modes [de- 
scribed by Eq. (j4.6b|) ]. their longest decay length is of 
order A^A = W . By choosing Lq = \QW we ensure that 
these modes are filtered out. 



filter 



lead 



filte 



L 



lead 



FIG. 4: Potential profile of a strip (length L), connected to 
leads by a pair of filters (length Lo). The Fermi wave vector 
in the leads is taken infinitely large; the finite Fermi wave 
vector in the filters removes the spurious evanescent modes 
excited by the leads. 



C. Conductance 

We have calculated the conductance at fixed Fermi en- 
ergy E = 2hv/ A as a function of the potential step height 
Uq. Results are shown in Fig. [5] for aspect ratio W/L = 3 
and lattice constant A = 10~^ L (solid curve) and com- 
pared with the solution of the Dirac equation (dashed 
curve). The agreement is excellent (for a twice smaller 
A the two curves would have been indistinguishable). 

The horizontal dotted line in Fig. [5] indicates the value 

lim lim (LlW)G/Go = I/tt (4.7) 

W/L^oo Uo-*E^ 

of the minimal conductivity at the Dirac point for a large 
aspect ratio of the strip. The oscillations which develop 
as one moves away from the Dirac point are Fabry-Perot 
resonances from multiple reflections at a; = and x = L. 
The filters of length Lq are not present in the contin- 
uum calculation (dashed curve), but the close agreement 
with the lattice calculation (solid curve) shows that the 
filters do not modify these resonances in any noticeable 
way. The filters do play an essential role in ensuring 
that the minimal conductivity reaches its proper value 
(|4.7p : Without the filters the lattice calculation would 
give a twice larger minimal conductivity, due to the con- 
tribution from the spurious evanescent modes of large 
transverse momentum. 



V. TRANSPORT THROUGH DISORDER 

We introduce disorder in the strip of length L by 
adding a random potential 5U to each lattice point, dis- 
tributed uniformly in the interval (— Ai7, AJ7). Since our 
discretization scheme conserves the symplectic symmetry 
exactly, there is no need now to choose a finite correlation 
length for t he p otential fluctuations (as in earlier numeri- 
cal studies [lll[il,[il,[i3,[ii,[2l[2l,[23). Instead we can 
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FIG. 5: Solid curve: conductance in the geometry of Fig. 2] 
The Fermi wave vector {E — Uo)/hv in the strip of length L 
and width = 3L is varied by varying the potential step 
height Uq at fixed Fermi energy E = 2?iu/A. The lattice 
constant A = 10~^ L. Dashed curve: the result from the 
Dirac equation (calculated from the formulas in Ref. [3^). 
corresponding to the limit A — > 0. The horizontal dotted line 
is the mimimal conductivity at the Dirac point. 
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FIG. 6: Same as Fig.U but now for the case that the potential 
in the strip fluctuates around the Dirac point: U = Uo + SU, 
with Uo = E and SU uniformly distributed in the interval 
{-AU,AU). 



let the potential of each lattice point fluctuate indepen- 
dently, as in the original Anderson model of localization 



A. Scaling of conductance at the Dirac point 

When Uq = E the potential Uo + 6U in the strip fluctu- 
ates around the Dirac point (see Fig. [6]). Results for the 
scaling of the average conductivity a = {L/W){G) with 
system size are shown for different disorder strengths in 
Fig. [T] We averaged over 3000 disorder realizations for 
L/A = 17,41,99 and over 300 realizations for L/A ~ 
239. The aspect ratio was fixed at W/L ^ 3. 

For sufficiently strong disorder strengths AU ^ 3hv/A 
the data follow the logarithmic scaling [ll|, [l2| 
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FIG. 7: Scaling with system size of the average conductivity 
a = {L/W){G} in a disordered strip at the Dirac point (ge- 
ometry of Fig. |6} . The length L of the strip is varied at fixed 
aspect ratio W/L = 3. The data are collected for different 
disorder strengths AU (listed in units of hv/A). 




a/Go = c\ii[L/l*{AU)] 



(5.1) 



10 100 1000 

L/r 



FIG. 8: Dependence of the conductivity of Fig. [7] on the 
rescaled system length L/l* (AU). The two dotted lines are 
the analytical weak and strong disorder limits. 



There is a consensus in the literature that c = I/tt can be 
calculated perturbatively [1^ as a weak antilocalization 
correction. The quantity I* plays the role of a mean free 
path, dependent on the disorder strength. We fit this 
scaling to our data with a common fitting parameter c 
(disregarding the data sets with low AU as being too 
close to the ballistic limit). The fitting gives I* for every 
data set with the same AU. 

The resulting single-parameter scaling is presented in 
Fig. [8] (including also the low AU sets, for complete- 
ness). The data sets collapse onto a single logarithmi- 
cally increasing conductivity with c « 0.33(1), close to 
the expected value c = I/tt « 0.318. To assess the im- 
portance of finite-size corrections [s^ we include a non- 
universal lattice-constant dependent term to the loga- 
rithmic scaling: a/Go = c\n[L/l*{AU)] + f{AU)A/L. 
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FIG. 9: Same as Fig. [S] but now for the variance of the con- 
ductance (instead of the ensemble average). The horizontal 
dotted line is the analytical prediction (15. 2|) . Var(G/Go) — 
0.116 W/L with W/L = 3. 



We then find c « 0.316(5), again close to the expected 
value [28] . These results for the absence of localization of 
Dirac fermions are consistent with earlier numerical cal- 
culations [m, [13] using a momentum space regularization 
of the Dirac equation. 



B. Conductance fluctuations at the Dirac point 

The sample-to-sample conductance fluctuations at the 
Dirac point were calculated numerically in Ref. [l6j us- 
ing the tight-binding model on a honeycomb lattice. An 
enhancement of the variance above the value for point 
scatterers was observed, and explained in Ref. @ in terms 
of the absence of intervalley scattering. A perturbative 
calculation of VarG = (G^) - (G)^ gives 



VarG 



3C(3) W 



Gl W/L:$>1. 



(5.2) 



Intervalley scattering would reduce the variance by a fac- 
tor of four, while trigonal warping without intervalley 
scattering would reduce the variance by a factor of two. 

In Fig. [9] we plot our results for the dependence of 
the variance of the conductance on the rescaled system 
size L/l*, with the AU dependence of /* obtained from 
the scaling analysis of the average conductance in Sec. 
IV Al The convergence towards the expected value (15. 2p 
is apparent. The numerical data of Fig. [9] supports the 
conclusion of Ref. [28] , that the statistics of the conduc- 
tance at the Dirac point can be obtained from metallic 
diffusive perturbation theory in the large-i limit. 

The tight-binding model calculation of Ref. [16] only 
reached about half the expected value (|5.2p , presumably 
because the potential was not quite smooth enough to 
avoid intervalley scattering. This illustrates the power of 
the finite difference method used here: We retain single- 
valley physics even when the correlation length of the 
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FIG. 10: Crossover from ballistic to diffusive conduction away 
from the Dirac point. The conductivity is plotted versus sys- 
tem size, at fixed Fermi wave vector {E — Uc)/Tiv = 0.8 
in the strip and fixed aspect ratio W/L — 3. The data is for 
difi^erent disorder strengths AU, listed in units of Tiv/ A. The 
dotted curves are a fit to the semiclassical formula (|5.3|l , with 
the transport mean free path Zo as a fit parameter. 




FIG. 11: Same as Fig. 1101 but now for the variance of the 
conductance. The data is plotted as a function of the rescaled 
sample size, using the values of the mean free path obtained 
from the fit of the conductance. The horizontal dotted line is 
the analytical prediction (|5.2|l . 



potential is equal to the lattice constant. 



C. Transport away from the Dirac point 



The results of Sees. IV Al and IV Bl are for potential 
fiuctuations around the Dirac point [Uq = E). In this 
subsection we consider the average conductance and the 
conductance fiuctuations away from the Dirac point. We 
take {E — Uq) — 0.8 hv/ A and vary the sample length 
L at fixed aspect ratio W/L = 3. The resulting size de- 
pendence of the conductivity is presented in Fig. \TU[ for 
different disorder strengths AU. 

Since antilocalization is a relatively small quantum 
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correction at these high Fermi energies, we are in the 
regime described by the semiclassical Boltzmann equa- 
tion [sl, [13] ■ In App. Owe apply a general theory ]4lj | 
for the crossover from ballistic to diffusive conduction, to 
arrive at the formula 



(G) - ^GoN^tri.J^^, 



(5.3) 



for the average conductance in terms of the transport 
mean free path Iq and the number A^strip = \E — 
UoKW/nhv) of propagating modes in the strip. From 
the fit of (G) versus L in Fig. [10] we extract the depen- 
dence on AU of lo, and then we use that information to 
investigate the scaling of the variance of the conductance 
with system size. As seen in Fig. [Til the variance scales 
well towards the expected value (|5.2[) . 



VI. CONCLUSION 

In conclusion, we have presented in this paper what 
one might call the "Anderson model for Dirac fermions" . 
Just as in the original Anderson tight-binding model of 
localization [s^, our model is a tight-binding model on 
a lattice with uncorrelated on-site disorder. Unlike the 
tight-binding model of graphene (with nearest neighbor 
hopping on a honeycomb lattice), our model preserves 
the symplectic symmetry of the Dirac equation — at the 
expense of a nonlocal finite difference approximation of 
the transfer matrix. 

Our finite difference method is based on a discretiza- 
tion scheme developed in the context of lattice gauge the- 
ory [26l [2^ . with the purpose of resolving the fermion 
doubling problem. We have adapted this scheme to in- 
clude the chiral symmetry breaking by a disorder po- 
tential, and have cast it in a current-conserving trans- 
fer matrix form suitable for the calculation of transport 
properties. 

To test the validity and efficiency of the model, we 
have calculated the average and the variance of the con- 
ductance and compared with earlier numerical and ana- 
lytical results. We recover the logarithmic increase of the 
average conductance at the Dirac point, found in numer- 
ical calculations that use a momentum space rather than 
a real space discretization of the Dirac equation [ll], [l2| . 
The coefficient that multiplies the logarithm is close to 
I/tt, in agreement with analytical expectations [2^. The 
variance of the conductance is enhanced by the absence 
of intervalley scattering, and we have been able to con- 
firm the scalingwith increasing system size towards the 
expected limit [l3| — something which had not been 
possible in earlier numerical calculations 16,] because in- 
tervalley scattering sets in before the large-system limit 
is reached. 

Our calculations support the expectation [2^ that the 
statistics of the conductance at the Dirac point scales 
towards that of a diffusive metal in the large-system limit. 
This would imply that the shot noise should scale towards 




FIG. 12: Scaling with system size of the Fano factor (average 
shot noise power divided by average current) in a disordered 
strip at the Dirac point (geometry of Fig. [6]). The length L 
of the strip is varied at fixed aspect ratio W/L = 3. The 
data are collected for different disorder strengths AU (listed 
in units of hv/A). The dotted horizontal line is the value 
F = 1/3 for a diffusive metal. The dotted curve is a fit to 
F — 1/3 -I- a[b + ln(I//r)]~^, included in order to indicate a 
possible scaling towards the expected value. 



a Fano factor F = 1/3 42J- Earlier numerical studies 
using the momentum space discretization found a 
saturation at the smaller value of _F = 0.295. Our own 
numerical results, shown in Fig. 1121 instead suggest a 
slow, logarithmic, increase towards the expected F = 
1 /3. More research on this particular quantity is required 
for a conclusive answer. 

We anticipate that the numerical method developed 
here will prove useful for the study of graphene with 
smooth disorder potentials (produced for example by 
remote charge fluctuations), since such potentials pro- 
duce little intervalley scattering. Intervalley scattering 
is absent by construction in the metallic surface states 
of topological insulators (such as BiSb [l^]). These sur- 
face states might be studied by starting from a three- 
dimensional tight-binding model, but we would expect 
a two-dimensional formulation as presented here to be 
more efficient. 
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APPENDIX A: CURRENT CONSERVING 
DISCRETIZATION OF THE CURRENT 
OPERATOR 

We seek a discretization of the current operator (|2.3[) 
that satisfies the condition (|2.6p of current conservation. 
Substitution of the expression (|2.19p into the condition 
(|2.6p gives the requirement 



(Al) 



The requirement that Eq. (jAl[l holds for any choice of po- 
tential fixes the discretization (|2.20p of the current oper- 
ator [up to a multiplicative constant, which follows from 
the continuum limit (|2.3p ]. 

This is an appropriate point to note that current con- 
servation could not have been achieved if the potential 
would have been discretized in a way that would have 
resulted in a nonsymmetric matrix Vm- For example, if 
instead of Eq. (|2.10p we would have chosen 

V'i' |(Kn+l,ri^'m+l,n + V"m+l,ri + l^m+l,n+l 

m.n + K^,„+l*m,„+l), (A2) 

with Vjn.n = V{xrmyn)i then the corresponding matrix 
Vm would have been asymmetric and no choice of Jx 
could have satisfied Eq. (|Aip . 



APPENDIX B: STABLE MULTIPLICATION OF 
TRANSFER MATRICES 



To perform the multiplication (|2.2ip of transfer matri- 
ces in a stable way (avoiding exponentially growing and 
decaying eigenvalues), we use the current conservation 
relation (|2.6p to convert the product into a composition 
of unitary matrices (involving only eigenvalues of unit 
absolute value). The same method was used in Refs. 
[ill [23 , H^, but for a different current operator, so the 
required transformation formulas need to be adapted. 

We separate the spinor degrees of freedom s = ± of 
the transfer matrix Aim into four N x N blocks, 



Mr 



M: 



M; 



(Bl) 



The current conservation relation (j2.6p with current op- 
erator (|2.20p can be written in the canonical form. 



Ml 



1 

-1 



M, 



1 

-1 



(B2) 



in terms of a matrix Mm related to Mm by a similarity 
transformation. 



Mr 



1/2 fJ'/^ J'/^ 



jl/2 _J^l/2 / • f^^^) 



Eq. (|B2]) follows only from Eqs. (HU and if the 

matrix TZ is Hermitian, which it is since is Hermitian 
with only positive eigenvalues [see Eq. (|2.14p ]. 



It now follows directly from Eq. (|2.6p that the matrix 
Um constructed from Mm by 

= (c d) ^ = (a - Id-' c bd-^) (B4) 

is a unitary matrix. Matrix multiplication of Mm^s in- 
duces a nonlinear composition of J7m's, 



M1M2 <^ C/i U2, 



defined by 



Al BA (A2 B2\ _ fAs B3 

Ci dJ^ {C2 D2J - {C3 D3 

As^ Ai+Bi{l-A2Di)-^A2C\, 

B3 = Bi{l-A2Di)-^B2, 

C3 ^ C2{1 - DiA2)-^Ci, 

D3 = D2 + C2{1- DiA2r^DiB2. 



(B5) 



(B6a) 

(B6b) 
(B6c) 
(B6d) 
(B6e) 



To evaluate the product (|2.2ip of MrnS in a stable 
way, we first write it in terms of the matrices Mm, 



M^n-nH Mm] n. 

\m=l / 



(B7) 



We then transform each transfer matrix Mm into a uni- 
tary matrix Um according to Eq. (|B4p and we compose 
the unitary matrices according to Eq. (|B6p . Each step in 
this calculation is numerically stable. 

At the end of the calculation, we may in principle 
transform back from the final unitary matrix U to the 
transfer matrix M = TiT^MTZ by means of the inverse 
of relation (lB4l) , 



[A B\ ,~, fC-DB-^A DB-^ 

^=[c d)^-^=[ -b-^a b-^ )■ 

This inverse transformation is itself unstable, but we may 
avoid it because [as we can see by comparing Eqs. (|B4p 
and (jBSp with Eq. (|3.1ip ] the final U is identical to the 
scattering matrix S between leads in the infinite wave 
vector limit. Hence the conductance can be directly ob- 
tained from U via the Landauer formula p.Sp [with the 
T„'s being the eigenvalues of BB^^ and CC^. 



APPENDIX C: CROSSOVER FROM BALLISTIC 
TO DIFFUSIVE CONDUCTION 

Away from the Dirac point (for Fermi wave vectors 
kp = \E — UqI/Tiv in the strip large compared to 1/L) 
conduction through the strip is via propagating rather 
than evanescent modes. If the number A^strip — kpW/ir 
of propagating modes is ^ 1, the semiclassical Boltz- 
mann equation can be used to calculate the conductance. 

As the transport mean free path Iq is reduced by adding 
disorder to the strip, the conduction crosses over from 
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the ballistic to the diffusive regime. How to describe this 
crossover is a well-known problem in the context of ra- 
diative transfer [iH . An exact solution of the Boltzmann 
equation does not provide a closed-form expression for 
the crossover, but the following formula has been found 
to be accurate within a few percent: 

(G) = CGoN^tripj^. (CI) 



The coefBcient Cd depends on the dimensionality d: C3 = 
4/3, C2 = 7r/2, Ci = 2. The length ^ is the socalled 
extrapolation length of radiative transfer theory, equal 
to Iq times a numerical coefficient that depends on the 
reflectivity of the interface at a; = and x = L. An 
infinite potential step in the Dirac equation has — Iq, 
see Ref. [42*1 . Substitution into Eq. (jCip then gives the 
formula (|5.3p used in the text. 
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